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Abstract 

We investigate a new numerical procedure to compute fermionic correlation func- 
tions at very small quark masses. Large statistical fluctuations, due to the presence 
of local "bumps" in the wave functions associated with the low-lying eigenmodes of 
the Dirac operator, are reduced by an exact low- mode averaging. To demonstrate 
the feasibility of the technique, we compute the two-point correlator of the left- 
handed vector current with Neuberger fermions in the quenched approximation, 
for lattices with a linear extent of L « 1.5 fm, a lattice spacing a ~ 0.09 fm, and 
quark masses down to the e-regime. By matching the results with the correspond- 
ing (quenched) chiral perturbation theory expressions, an estimate of (quenched) 
low-energy constants can be obtained. We find agreement between the quenched 
values of F extrapolated from the p-regime and extracted in the e-regime. 
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1 Introduction 



One of the most important challenges for lattice QCD is the simulation of quarks 
with masses light enough to reach kinematical regions where chiral perturbation the- 
ory (ChPT) ^ [2j can be verified and safely applied. Despite valuable efforts in the 
infinite- volume limit [31-0) it is still unclear how small the quark masses need to be 
in order to reach the chiral region (for recent developments, see [7||S]). More recently 
the e-regime of QCD ^ in which m — > at finite volume V, is being attacked 
on the lattice El El 1^1 El El E] ■ Difficulties to reach these kinematical corners 
arise because most of the numerical techniques currently used on the lattice become 
inefficient when the quark masses approach the chiral limit. 

In the presence of spontaneous chiral symmetry breaking, chiral perturbation the- 
ory suggests that the low-lying eigenvalues of the massless QCD Dirac operator scale 
proportionally to {T,V)~^, where S is the bare quark condensate and V is the lat- 
tice volume "RIEI- In addition random matrix theory |18| I19j (and to some extent 
quenched ChPT [201) predicts the distribution of each individual eigenvalue as a func- 
tion of T,V. Comparisons of these predictions with numerical data obtained on the 
lattice in quenched QCD have been reported by several collaborations |211 llll [TH I15j . 
Only recently has an extensive study of this issue, performed at several volumes and 
lattice spacings, shown a detailed agreement between some of these predictions and 
quenched QCD for volumes larger than about 5 fm^ jl5j . No predictions are available 
for the properties of the corresponding wave functions (for a recent compilation of nu- 
merical results, see [HI). It is an empirical observation, though, that they can develop 
local "bumps" with a non-negligible probability ^^1 ■ 

It is then conceivable to exploit such basic properties of QCD as suggested by 
ChPT to develop exact and generic numerical algorithms which remain efficient when 
the quark masses approach the chiral regime. For instance, due to the scaling of the 
low-lying eigenvalues with V, the computation of the fermion propagator with stan- 
dard techniques becomes very demanding when the quark masses approach the chiral 
limit, since the Dirac operator tends to be strongly ill-conditioned. The slow-down can 
be dramatically reduced by subtracting a few of the lowest modes and treating them 
separately f23] . 

In this paper we propose a new technique^ to compute QCD fermionic correlation 
functions at very small quark masses on the lattice. An exact low-mode averaging 
procedure is introduced to reduce large statistical fluctuations induced by the pres- 
ence of "bumps" in the wave functions associated with the low-lying eigenvalues of the 
Dirac operator which happen to occur where the fermionic operators are localized. The 
feasibility of the technique is proven by computing the correlator of two left-handed 
vector currents in quenched QCD with Neuberger fermions |24j-|3nj. We find that the 
variance of the estimate is markedly reduced with respect to the one with standard 
techniques when the quark mass is decreased, and the e-regime can be safely reached 
in all topological sectors. 

The use of the quenched approximation mainly serves to test our method, which 

^During the writing of this paper, T. DeGrand and S. Schaefer applied a technique similar to the one 
presented in this work to reduce statistical noise in the computation of two-point correlation functions 
of bilinear operators |31|. A similar idea has also been sketched by R.G. Edwards to reduce noise in 
the computation of singlet correlation functions \'22]. 
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we expect will be effective for simulations of full QCD as well. It is known that the 
quenched theory is afflicted with several problems. In particular, the removal of the 
fermion determinant renders the theory non-unitary. An effective low-energy descrip- 
tion of the quenched theory is formally obtained if an additional expansion in 1/Nc, 
where N^. is the number of colors, is carried out together with the usual one in quark 
masses and momenta. The resulting so-called quenched ChPT [321 133j leads to infrared 
divergences in certain correlation functions. These divergences reflect, at least par- 
tially, the sickness of quenched QCD. Here we adopt the pragmatic assumption that 
- despite the fact that it is not an asymptotic expansion of quenched QCD (at fixed 
Nc) - quenched ChPT describes the low-energy regime of quenched QCD in certain 
ranges of kinematical scales, where correlation functions can be parameterized in terms 
of effective coupling constants, the latter being defined as the couplings which appear in 
the Lagrangian of the effective theory. With this assumption in mind we compare the 
predictions of quenched ChPT with the numerical results obtained in the kinematically 
accessible ranges in the p- and e-regimes for the correlator of two left-handed currents. 
Thereby we estimate values of F and 05 (to be defined below), which, according to our 
working assumption, are then identified with quenched versions of the physical LECs. 

It is interesting to note that the correlator of two left-handed vector currents is free 
from zero-mode contributions and, at fixed volume, remains finite when the quark mass 
m — > 0. This is in contrast to the correlators of two scalar or pseudoscalar densities in 
sectors of non-zero topological charge at finite volume. They develop divergences 
with residues given by correlation functions of zero- mode wave functions |16j . Since in 
both cases the leading non-trivial behaviour is governed by F, these two different types 
of correlators offer independent determinations of this coupling. 

The paper is organized as follows: in Sect.Elwe collect the formulae for the two-point 
correlation function of the left-handed vector current in chiral perturbation theory, in 
Sect. 01 we describe the low- mode averaging procedure for fermionic correlation func- 
tions, and in Sect. we report the details of the simulations we have performed in 
quenched QCD, as well as the results obtained. We conclude in Sect. El In the first 
Appendix we provide more details of our notations, and in the second we collect further 
useful formulae obtained in ChPT at finite volume. 

2 Left-handed current correlator in ChPT 

We start by considering the physical, unquenched QCD. The Euclidean Lagrangian 
of the chiral effective theory for this case reads, at leading order in the momentum 
expansion, 

C = ^Tt {d^U^d^u} - I Tr {e^^/^ff/M + Aft^/te^^/^f } , (1) 

where U G SU(A^f) is the meson field, A^f = 3, Af = diag(m„, m^, m^) the quark mass 
matrix, 9 the vacuum angle, and, in the chiral limit, F equals the pseudoscalar decay 
constant and S the chiral condensate At the next-to-leading order (NLO) in the 
momentum expansion, additional operators appear in the effective Lagrangian, with 

^In this section we use continuum notation. 
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the associated LECs ai,a2,--- H] ^. In the fohowing we consider only the case of a 
degenerate mass matrix, i.e. niu = ma = mg = m. 

The LECs can be determined by computing suitable QCD correlation functions 
on the lattice at small masses and momenta, and by comparing the results with the 
parameterization given by ChPT. In this paper we are interested in the two-point 
correlation function |231 134j 



of the left-handed charge density Jq{x), where t 
leading order in momentum expansion, it reads 

F 



d'x{js{x)j^m (2) 

xq. In the effective theory, at the 



(3) 



where T"^ is a traceless generator of SU(A'f), acting on flavor indices, and on the side 
of QCD we assumed the normalization of Eq. H20|l below. 

A kinematical range of scales which is suitable for extracting the LECs in a fi- 
nite box of volume V = TL^, is the so-called p-regime. It is defined by the con- 
straints MpL > 1 and Mp <^ 4ttF, where Mp is the pseudoscalar meson mass. Writing 



C^it) = Tr [T"T'']C(i), 
be expressed as jSH] 



the next-to-leading order (NLO) finite-volume prediction can 
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and = 2mT,/F'^. The effective finite- volume meson mass Mp and the functions gi 
and G{M'^) are reported in Appendix B. Finite-volume effects are exponentially small 
if MpL ^ 1, and we can set gi = 0. 

A less explored kinematical region of QCD, where the LECs can be extracted, is 
the so-called e-regime, where MpL <C 1 and the linear extent of the box is such that 
AnFL 3> 1 [01 ^Ij • In this regime topology plays an important role , and for fixed 
topological charge v jSHl EZl OS 
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^We adopt the convention of |35( where = 8(47r)^Li, with Li as defined in 
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The constants /?i and /cqo are related to the (dimensionally regularized) value of 

by 

Furthermore crj^(/i) = A''f^^d{lndet[/,y+j_i(//)]}/d//, where the determinant is taken 
over an A'f x A'^f matrix, whose matrix element is the modified Bessel function 

Iu+J-^ isHiini. 

In ChPT there is a well-defined prescription to compute correlation functions at 
fixed topology, and here we assume that also in QCD they have a well-defined meaning 
in the continuum limit at non-zero physical distances. Although plausible, this is a non- 
trivial dynamical issue and to pose precise questions we must introduce an ultraviolet 
regularization. We here adopt a lattice regularization with fermions discretized d la 
Neuberger [^. The massless Dirac operator obeys the Ginsparg- Wilson (GW) relation 
and therefore preserves an exact chiral symmetry. The topological index assigned to a 
configuration is = n+ — n_, where n_|_ (n_) are the numbers of zero- modes of D with 
positive (negative) chirality. Our working hypothesis is that correlators of composite 
operators at non-zero physical distance have a continuum limit in any given sector of 
index u independent of the particular choice of ^. Some recent numerical evidence 
(in the quenched approximation) consistent with this scenario can be found, e.g. in 

Refs. [iniinn]. 

In quenched ChPT the flavor singlet field cannot be integrated out and therefore 
additional coupling constants appear in the chiral Lagrangian |321 133j . In particular 
the singlet mass parameter m'^/2Nc is dimensionful and, even if suppressed by large- 
Nc counting, cannot be tuned by changing the kinematical conditions. Consequently, 
the standard chiral expansion is expected to be useful only in a window of Euclidean 
momenta , where vn^jlNc <C (^^ <C (47rF)^. 

In the quenched theory, Eq. Q remains the same, apart from the omission of the 
last constant term oc A'^f, but at NLO the interpretation of the parameters in terms of 
those of quenched ChPT changes flQ] ^ . In particular, the parameter Fp is volume- 
independent, 

= Fp = F 

where the LEC 05 is finite at this order. In the e-regime, on the other hand, the 
correlator of two left-handed vector charges is modified to be |37l I34j 

CW = ^{l + ^-.(M)/^i(^)}, (12) 

where in this case j41| 

cj^(/i) = /i{/^(/x)K^(/i) + 4+i(/i)AV_i(/i)} + ^ , (13) 

''Since the space of lattice gauge fields is connected, different choices of D possibly lead to different 
assignments of index for a given configuration. 

^Analogous LECs in ChPT and quenched ChPT are indicated with the same symbols, since they 
can be clearly distinguished from the context. 
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and ly and K^, are modified Bessel functions. 

3 Low-mode averaging 

Although the technique we are going to describe can be more widely applied, we re- 
strict ourselves to lattices of spacing a, volume V = TL^ and with periodic boundary 
conditions imposed on all fields. QCD gluons and fermions are discretized using the 
standard plaquette Wilson action and the Neuberger-Dirac operator D, respectively. 
The latter pS] satisfies the Ginsparg- Wilson relation [25| 

75£» + D75 = a^75^ , (14) 

and thus preserves an exact chiral symmetry at finite lattice spacing _29j. The Neu- 
berger operator, the parameter a and other conventions used in this section are defined 
in Appendix A, and are the same as in Ref. The massive lattice Dirac operator is 
given by 

CiTYl 

Dm = {l-—)D + m (15) 

where < am < 2. For a given gauge configuration the massless Dirac operator can be 
diagonalized, and chirality implies that non-chiral modes appear in complex conjugate 
pairs, i.e. 

Drjx^ = AfcT/Afc, k = l,2,..., (16) 
Dv\l = Kmi , V\l = 75??Afe • (17) 

Random matrix theory jl8l I19j (and to some extent quenched ChPT [20] ) predicts the 
probability distribution of each eigenvalue in the low-lying end of the spectrum to be a 
function of the rescaled variable C = | A|Sy only. Some of these predictions are reported 
in Fig. ^ for the quenched theory. In particular, the predicted individual distributions 
of the first four eigenvalues together with the total microscopic density for = are 
shown in the plot on the left, while the individual distributions of the first non-zero 
eigenvalue are shown for several values of v in the right plot. It is interesting to note 
that there is no gap for small values of (, and the total distribution is predicted to be 
Ps{C) ~ ^2(|!/|+Aff)+i^ Therefore arbitrarily small eigenvalues can occur (either in the 
full or the quenched theory) with a probability which decreases exponentially with 
and A^f. The expectation value of the lowest eigenvalue and the level splittings near 
the origin are of 0{1/T,V), as can be seen in Fig. ^ The ratios of expectation values of 
low-lying eigenvalues are then parameter-free predictions of RMT, and they have been 
confronted with quenched QCD data in Ref. ^Sj. An impressive agreement has been 
found for volumes larger than about 5 fm^. 

No analytic predictions are available so far for the properties of the corresponding 
wave functions. In Ref. jl5j it was found that the probability distribution of their norm 
at a fixed lattice site is broader than the one expected for a normalized random vector^ . 

®The behaviour of the eigenfunctions corresponding to the low modes of the Dirac operator has 
already been studied numerically in different contexts, see Ref. |22| . 
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Figure 1: Probability density distributions of low-lying eigenvalues in quenched QCD. 
Left: distributions of the first four eigenvalues Pi(C) {i = 1; • • • 4) and of the microscopic 
total density Ps{C) for = 0. Right: distributions of the first eigenvalue for =0, 1, 2. 



Since in the following we are interested in the correlation function of two left-handed 
currents, we have studied the probability distribution of u = |Re[t>ii(x)]|y, where 

Vkiix) = ivljoP-rixMx) , (18) 

on a lattice of volume 16^ at /3 = 6.0 (Lattice B, see next section). The result is shown 
in Fig. [21 together with the analogous prediction for the case that the r/'s are treated 
as random vectors with unit norm and such that Y^x[''lx^l5'n\k]i^) — can be 

seen from the inset in the figure, the distribution decreases quite slowly, and thus the 
probability of finding a point on the lattice for which u exceeds its mean value by far is 
quite high: for instance, the probability for finding u > 0.5 is 0.7%. We also mention 
that the distribution of |Im[t;ii(x)]|y nicely overlaps with the one reported in Fig. |2j 

These properties of eigenvalues and eigenfunctions of the Dirac operator imply 
that, when a fermionic correlation function is computed with standard Monte Carlo 
techniques, the relative contributions from the various eigenvalues can change dramat- 
ically with m: 

• for m ^ IfEV, the low-lying end of the spectrum of Dm is very dense near m, 
and (many) contributions from the corresponding wave-functions are averaged 
with essentially the same weight. 

• for m ~ l/TiV, the mass is comparable with the expectation value of the modulus 
of the lowest eigenvalue of D (see Fig. and therefore the low-lying spectrum of 
Dm appears discrete (the splittings of the eigenvalues are of the same order as their 
values). The contribution of just a few eigenvectors to a given observable can be 
substantial and their space-time fluctuations can then induce large fluctuations in 
its estimate. A significant improvement can be achieved in this case by low-mode 
averaging as we explain in the remainder of this section. 
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Figure 2: Probability distribution of n = |Re[uii(x)]|y (cf. Eq. (|18j) ) for the first norma- 
lized non-zero mode eigenvector of D. The solid line is the distribution obtained at 
(5 = 6.0 on a 16^ lattice (Lattice B, see below), while the dashed curve is the one 
expected for a random vector model. 



• for m <C 1/TjV , the mass is much smaller than the expectation value of the lowest 
eigenvalue of D. Extremely small eigenvalues of Dm can then occur with a small 
but non-negligible probability. The standard Monte Carlo sampling of the path 
integral for fermionic correlation functions is problematic in this case, because a 
sizeable contribution to the estimate of these integrals is obtained from configura- 
tions that have a small statistical weight (i.e., those with the smallest eigenvalues 
of Dra, for which these observables are very large). We did not find a straightfor- 
ward solution to this problem. It would probably require an improved algorithm 
where the low-mode contribution to the observable is included in the Monte Carlo 
sampling of the integral. If such an algorithm exists, it will presumably be very 
expensive in practice. 

The information above can be exploited in order to develop efficient algorithms for 
computing fermionic correlation functions in the regime m>l/TiV. In the following 
we focus on the two-point function 

c^\t) = Y,mMm (19) 

X 

of two left-handed charge densities, 

J^{x) = il^{x)TV-i'{x) , (20) 




7 



where -0 is defined in Appendix A. Writing C"'*(t) = Tr [T''T*]C(t), and using the 
spectral decomposition, we get 

= - E ( E ''^ [vlloP-Vx^{x)[viloP-VxM (21) 



k,l 



EE 



where = {(1 — am/2)Xk + m} are the eigenvalues of the massive operator Dm- 

The standard way of estimating C{t), by computing for every gauge configuration 
the propagator from a local source to any other point, turns out to be efficient in 
the mass range m S> 1/TiV, where the correlator is the result of an average over 
many comparable contributions from different eigenfunctions. When the quark mass 
reaches the e-regime, i.e. m ~ IfEV, the presence of "bumps" in single wave functions, 
associated with low-lying eigenvalues, can generate much larger statistical fluctuations. 
They can be reduced by noting that each contribution to the sum over k and / in Eq. H22|) 
can be estimated independently. The statistical error of each term can be reduced by 
increasing the number of configurations, exploiting the symmetries of the theory, etc. 

In particular, the variance of the low-mode contribution to C{t) can be decreased by 
exploiting the translational invariance of the theory. It is important to note that, based 
on the techniques developed in Ref. j23j . this can be achieved without computing low- 
lying eigenmodes and eigenfunctions with very high precision, which can be numerically 
expensive, especially with Neuberger fermions. 

We now describe a particular implementation of the ideas sketched above, which we 
have adopted for computing C (t) . For each gauge configuration we have computed the 
topology v, and, starting from a set of Gaussian random sources, we have minimized 
the Ritz functional until the estimated relative errors of the calculated eigenvalues drop 
below a specified bound ujk 123], cf. Eq. (|23j) . This guarantees that the orthonormal Dirac 
fields ni • • • m„, which approximate the eigenvectors, satisfy 

Auk = akUk + rk, {ui,rk) = 0, \\rk\\ < ujkOk , (23) 

for all k, I. In Eq. (|23() . A = PsDmDmPs, Ps {s = ±) is the projector into the chirality 
sector without zero modes, and are the approximate eigenvalues of A. We can 
then define a "subtracted" left-left propagator P-S^{x, y)P+ for the massive Neuberger 
operator (which we have computed in practice as described in Ref. j23j ) through the 
equation 

P.S{x,y)P+ = P\y2—ek{x)ekiy)^ +S''ix,y)]p+ (24) 



where 



aui 



d^ 



S{x, y) = ^ _ ,J D-'}ix, y) (25) 
1 — am/ 2 



Ck = PsUk + P^sDPsUk . (26) 



'^Unlike tja^ , the two chiral components of are not normalized, in order to simplify our formulae. 
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The two-point correlation function in Eq. ()19() is then given by 

C{t) = C7" (t) + C^' (t) + C'^^ (t) , (27) 



where 



k,l=l X 



n 

C"'(t) = el(x)7oP-5^x,0)7oP-efc(0)) + (x^O), (29) 



fc=i 



= -Y^{TT[-foP-S\x,OhoP-SHo,x)]) , (30) 

X 

and the space-time dependences of eigenvectors and propagators have been shown ex- 
phcitly. By noticing that the starting vectors of the Ritz minimization procedure are 
extracted with a translationally invariant action, it is straightforward to prove that 
each contribution on the right-hand side of Eqs. (|28j) - (|,''{()|l is translationally invariant 
even if the vectors in Eq. (|23|) are only approximate eigenvectors of A, i.e. ^ 0. 
In this case, in addition to the gluon field, also the random vectors needed to start the 
Ritz functional minimization should be translated^. Therefore, 



k,l=l x,y 



^"'W = 4(^)7o^'-5'^(x,y)7oP~efc(y)) + (x^y) (32) 

k=l x,y 

hold independently of the number n of eigenvectors which have been subtracted and of 
the precision oj^ they have been determined with. By contrast, the statistical variance 
of the signal changes with n and LOk- Note that the computation of C^\t) can be quite 
expensive since it requires an inversion of the Dirac operator for every vector e^. 



4 Numerical results 

To test the procedure described in the previous section, we have simulated two lattices 
at P = 6.0 (a ~ 0.09 fm) with volumes V = 24 x 16^ (A) and V = 16^ (B). Since the 
generation of gauge-field configurations consumes a negligible amount of computer time, 
we have performed many update cycles between subsequent measurements so that they 
can be assumed to be statistically independent. The computation of the index, the low- 
lying eigenvalues and the inversion of the Neuberger-Dirac operator have been carried 
out using the techniques reported in Ref. 23 . For each gauge-field configuration of 
the lattice A(B), 7(6) low-lying eigenvalues of D^D have been extracted in the chirality 
sector without zero modes with a relative uncertainty of uJk = 0.05 (see below), while 
a further one has been determined with lower precision to stabilize the Ritz functional 

*We thank M. Liischer for having clarified this point to us. 
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minimization. The subtracted and the fuh propagator have been computed by requiring 
a residue of 5 • 10~^ in the adaptive conjugate gradient. 

Lattice A has been devoted to studying the left-left correlation function in the 
p-regime: we have generated 113 gauge configurations, inverted Dm for masses am = 
0.025, 0.040, 0.060, 0.080, 0.100, and computed the correlation functions in the standard 
manner (local source and sink) and with the low-mode averaging (LMA) procedure 
described in the previous section. After symmetrizing the correlators around t = r/2, 
we estimated the statistical errors by a jackknife procedure and fitted the correlation 
functions with the expression given in Eq. (jU in the time interval 6—12. The lower limit 
was fixed at the point where we found stabilization of the effective meson masses. The 
results of the fits are given in Tabled They are compatible with previous computations 
[421 I43j (for recent reviews see |441 I45| ') within the statistical uncertainties. In this 
regime and at this volume, the benefit of the low- mode averaging is visible for the two 
lightest quark masses only, and it is likely to be less effective when the volume increases 
and the quark mass is kept fixed. 

The values of Fp follow a remarkable linear behaviour in the quark mass. A linear 
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aF^P 
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aF^P 


0.025 


0.199(6) 


0.0341(6) 


0.198(8) 


0.0336(9) 


0.040 


0.242(5) 


0.0355(6) 


0.244(7) 


0.0349(7) 


0.060 


0.292(5) 


0.0374(6) 


0.295(6) 


0.0369(6) 


0.080 


0.335(4) 


0.0392(6) 


0.339(5) 


0.0390(5) 


0.100 


0.375(4) 


0.0410(6) 


0.380(5) 


0.0410(5) 



Table 1: Meson masses and decay constants from the left-left correlators computed in 
the p-regime with the LMA procedure, and with the standard local one. By fixing 
the lattice spacing with rg from Ref. [IS]) the physical value of the Kaon mass Mk = 
496 MeV, expressed in lattice units, would correspond to aMx ~ 0.234. 

fit of the form aFp = Ai + A2 ■ {am) gives 

Ai = 0.0318(7) , A2 = 0.093(6) LMA (33) 
Ai = 0.0310(10) , A2 = 0.100(10) Local . (34) 

Quadratic fits to the data give results very well compatible with the previous ones, with 
the coefficients of the quadratic terms compatible with zero. 

Lattice B has been reserved for the e-regime: we have generated 203 independent 
gauge configurations of which 31, 66 and 44 have topological charge = 0, 1, 2, re- 
spectively. In these topological sectors we have computed the quark propagators for 
masses am = 0.001,0.005,0.010 which correspond to mTiV ~ 0.11, 0.55, 1.1 respec- 
tively, if the bare quark condensate S is taken from the analogous lattice of Ref. |15j . 
As in the previous case, the two-point correlators of the left-handed current have been 
computed in the standard manner and with LMA. 

In Fig.|31the Monte Carlo histories for the left-left correlator at t/a = 6 and for the 
lowest non-zero eigenvalue of D^D are reported for \v\ = 0, 1 and am = 0.005, 0.001. 
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Figure 3: Monte Carlo histories (Lattice B) for the absolute value squared of the lo"west 
non-zero eigenvalue of D (bottom) and for the left-left correlators computed at t/a = 6: 
= (left) and \u\ = 1 (right), am = 0.001 (middle) and am = 0.005 (top). The 
dashed lines in the eigenvalue plots indicate the two values of (am)^. In the plots of 
the correlation functions, the thick lines are obtained "with LMA, while the thin ones 
in the standard manner. 



For the lightest mass, a spike in C{t) is clearly visible in correspondence with a very low 
(the lowest produced) eigenvalue, which happens to be roughly one order of magnitude 
lower than its expectation value |15j . A closer look at this configuration reveals that the 
spiky contribution is indeed due to the light-light contribution C^\t) and is not cured 
with the LMA procedure proposed here. As expected (cf. Sect.|2I), for m <C 1/TjV the 
Monte Carlo history shows evidence for extreme statistical fluctuations, and therefore 
we discard data at the lightest mass in the following analysis. The spiky behaviour 
disappears for the two heavier masses and for them the Monte Carlo history is well- 
behaved. 

For masses m ~ l/SV, the LMA estimate of the correlation function C{t) is indeed 
less fluctuating than the one computed in the standard manner. Its variance turns out 
to be roughly a factor two smaller for the topologies and masses that we consider, as 
shown in Fig. |l]for \v\ = 0, 1 and am = 0.005. The variance reduction is a function of 
the number of eigenvalues extracted n, and the precision oo^ they have been computed 
with. It can also vary depending on which contributions are chosen to be averaged 
over. We tried several values of n and lo^ and the ones used in this computation turn 
out to be a good compromise between the gain in statistics and the additional compu- 
tational cost (roughly a factor 2), which is mainly due to the computation of the mixed 
contribution C^\t). A more systematic optimization of these parameters is desirable 
but goes beyond the scope of this paper. 



11 




Figure 4: Left-left correlators in the e-regime for am = 0.005, |z^| = (left) and 
|z^| = 1 (right), computed with LMA (filled diamonds) and in the standard manner 
(open squares). The dashed lines represent fits to LMA data (see main text). 



We fitted the LMA correlations with the expression given in Eq. (|12|1 . 



C{t) = ^ 



1 + ^2 /ii ( 4 



(35) 



for each topological sector and mass. The lower temporal limit was fixed at t/a = 5, 
a point where we found stabilization of the for all correlators. The results obtained 
and the associated errors computed with a jackknife procedure are reported in Table |21 
Within the large statistical errors, the values for B2 are compatible with the prediction 



am 




Bi 


B2 







0.033(3) 


3(2) 


0.010 


1 


0.032(2) 


5(1) 




2 


0.030(2) 


5(1) 







0.034(4) 


2(2) 


0.005 


1 


0.032(2) 


4(1) 




2 


0.030(2) 


4(1) 



Table 2: Results from the e-regime data of a fit of the form given in Eq. (|35() . 



of Eq. ()12[) and the estimate of the condensate extracted in Ref. jl5j . The results for 
Bi are very well compatible for all masses and topological sectors. We estimate our 
best value of the decay constant in the e-regime, aF = 0.0312(11), by averaging the two 
values of Bi in each topological sector and then averaging the results as independent 
determinations. This value is in very good agreement with the linearly extrapolated 
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0.048 



0.044 



0.04 



0.036 



0.032 



0.028 



0.02 



0.04 



0.06 



0.08 



0.1 



0.12 



Figure 5: Results obtained in the p-regime for Fp (open circles) and in the e-regime 
for F (filled diamond). The dashed line represents a linear fit to the p-regime points. 
The corresponding intercept is also shown (filled circle). 



value from the p-regime, as shown in Fig. El and constitutes one of the main results of 
this paper. 

By renormalizing the local left-handed current with Zj = 1.55 j42j and by fixing the 
lattice spacing with ro from Ref. gHl, we obtain F = 104(2) MeV and F = 102(4) MeV 
in the p and e-regimes, respectively. These values are compatible within 2a yet more 
precise than the one we found recently in Ref. ^H] from a study of the topological zero- 
mode wave functions. Note however that the systematic error due to the finite lattice 
spacing has not been quantified in the present study. The agreement between these 
two determinations would provide a further check of our main working assumption 
that quenched ChPT describes the low-energy regime of quenched QCD in certain 
ranges of kinematical scales. By fitting Fp linearly as a function of (Mp)^, we also 
obtain our estimate of the LEG in Eq. (|1H) to be as = 1.66(8). This value is in the 
same ballpark as the one found in Ref. while it is higher than the one obtained 
in Ref. [SSI- '^o understand the discrepancy we have analyzed our data as suggested 
in Ref. j35j . and we have found as = 1.08(5), which agrees well with their value of 
0.99(6). Closer inspection shows that the difference can be traced back to two features 
of the analysis presented in |SSj: the first is the Taylor expansion of 1/(1 + as2/ref/2) to 
leading order in as. Although such an expansion is justified, since the difference with 
the exact result is a higher-order effect, the resulting ambiguity is large and produces 
an estimate for as which is smaller by 20-25%. Indeed, this discrepancy was included 
as a systematic uncertainty in j^. The remaining difference can be explained by the 
fact that in Ref. [HS] the value F was constrained to coincide with the experimental 
value of the pion decay constant, which is smaller by about 15% than typical quenched 
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estimates. This assumption is not in accord with our working hypothesis that quenched 
data should be reproduced by quenched ChPT. 

The errors quoted above for the LECs include only our statistical errors. A more 
detailed assessment of the various systematic errors would require computations at 
different volumes and lattice spacings, which goes beyond the scope and primary goal 
of this study. 

5 Conclusions 

The low-mode averaging technique proposed in this paper reduces large statistical fluc- 
tuations in correlation functions due to the presence of local "bumps" in the wave 
functions associated with the low-lying eigenmodes of the Dirac operator. When ap- 
plied to the two-point function of the left-handed vector current in the region of quark 
masses m ~ 1/Y,V, it provides an estimate of the correlator with a variance signifi- 
cantly reduced with respect to the standard one. As a result the e-regime of QCD can 
be safely reached in all topological sectors. 

It is conceivable that more involved correlation functions such as singlet diagrams 
(see for example |49| ) or those needed for non-leptonic weak decays can also benefit 
from low-mode averaging. Conventional formulations of lattice QCD may profit as well 
from this technique. For instance, it could mitigate the problem of exceptional config- 
urations encountered for Wilson fermions, or speed up unquenched simulations. 

By matching the quenched QCD results for the left-left correlators computed for 
quark masses in the p- and e-regimes with those of quenched ChPT, we obtained es- 
timates for the quenched low-energy constants F and a^. The agreement we found 
between the value of the pseudoscalar decay constant extrapolated from the p-regime 
and the one extracted directly in the e-regime is remarkable. 

In Ref. [2]) we studied the possibility of using the contribution of topological zero- 
mode wave functions to the pseudoscalar correlator, in order to extract F. As far as 
the convergence of the chiral expansion is concerned, it was found that in full QCD 
(assuming F ~ 93 MeV), one would need to go to a lattice extent > 2.0 fm (> 2.5 fm) 
to have the first non-trivial relative correction to be less than 50% (30%). Inspecting 
the constant part of the e-regime expression in Eq. ((7| , we find that for the observable 
studied in the present paper, the same relative corrections would be obtained already 
with lattice extents > 1.4 fm (> 1.8 fm). Therefore, current correlators near the chiral 
limit should allow for smaller systematic uncertainties in the extraction of the LEC 
F at realistically accessible volumes than the zero- mode wave functions. It would be 
interesting to study whether the same remains true for other LECs as well, such as 
those related to weak decays. 

In the quenched theory, on the other hand, systematic errors are very hard to quan- 
tify, but the apparent convergence of our expression for the current correlator compares 
well with the one observed for the zero- mode wave functions in Ref. J^] , while having at 
the same time the advantage that the singlet parameters ttiq/INc, a/2Nc of quenched 
chiral perturbation theory do not enter at all at this order. 
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Appendix A. Some definitions for the GW fermions 

In this paper we employ the same conventions as in Ref. [ISl- The Dirac matrices 
satisfy 

(7m)^ = 7m, {ifiHu} = 26f,^, (36) 
and we have chosen a chiral representation with 

75 = 70717273 = f Q j • (37) 



The chiral projectors are defined as 



P± = ^(l±75). (38) 



The Wilson-Dirac operator is given by 



= ^ {7m(V* + V^) - aV* V^} , (39) 



where 



Vf,ij{x) = -{U{x,fi)'ilj{x + afi)-ij{x)} , (40) 
a 

V;i;{x) = ^{i;{x)-U{x-afi,fx)~'i;{x-afi)} (41) 

are the gauge-covariant forward and backward difference operators, a denotes the lattice 
spacing, U (x, fi) G SU(3) are the link variables and fi is the unit vector along the 
direction /i. The Neuberger-Dirac operator is defined as 

D = i{l + 75sign(Q)}, (42) 
a 

Q = 75 (aZ^w - 1 - s) , a=-^, (43) 

1 + s 
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where s is a real parameter in the range |s| < 1. It satisfies the Ginsparg- Wilson 
relation in Eq. (|14() . Infinitesimal chiral transformations of the fermion field are given 
by [221 

'^V' = 75(1 - ! ^-f/i = ■075 . (44) 

The modified fermion field 

^ = (1 - \aD)i^ (45) 

transforms according to 

,5^ = 75-(/i , (46) 

and therefore if a composite operator is defined using ip instead of ij:, it has the same 
transformation behaviour as the corresponding one in the continuum. 



Appendix B. Pseudoscalar mass in the p-regime of ChPT 
to NLO 

For completeness we report in this appendix NLO expressions for the pseudoscalar 
meson mass in the p-regime of ChPT. We again consider the case of degenerate light 
quarks only. The effective finite- volume pion mass Mp , entering the prediction for the 
correlation function of two left currents reported in Eq. is given by 



Mj, EE Mp 1 + 



91 



where at the same order the infinite-volume mass Mp is 



1 + 



Mp = M 

The function gi reads |1H1 
1 dX ^ 



G{M 



2A^fF2 2(47rF)2 



N{a4, + 05 — 2NfaQ — 2as 



91 



(4vr)2 7o A2 



El 



and, in dimensional regularization ^, 



^nfi ) exp 



1 

'4AV 



(47) 



(48) 



(49) 



(50) 



In the quenched approximation (to the extent that it is well defined for this observable) , 
these predictions are modified to be 



Mp 



1 + 



1 



2Fl \ 2N, 



a ^ aMp — rriQ d 



2Nr 



dMl 



91 



with the infinite-volume mass Mp given by 



Mp = M 



1 + 



1 



2^2 \2Nc 



-G{M') + 



2Nc 



-H{M^ 



M' 



2(47rF)^ 



(a5 - 2as) 



(51) 



• (52) 



^The divergence of G{M^) for d ~ 4 cancels against those in the a^'s 
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In Eqs. 1)511) and (|52() . ajlN^ and rn^/2Nc are the parameters related to the flavour 
singlet field (with the normalisation conventions of [341 13 7j ) , and 



H{M'') = I ^\ , ^ ^ . (53) 
^ ' ' (27r)'^ (^j2 + M2)2 ^ ' 
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